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Abstract. 

Atoms, molecules or excitonic quasiparticles, for which excitations are induced by 
external radiation fields and energy is dissipated through radiative decay, arc examples 
of driven open quantum systems. We explain the use of commutator-free exponential 
time-propagators for the numerical solution of the associated Schrodinger or master 
equations with a time-dependent Hamilton operator. These time-propagators are 
based on the Magnus series but avoid the computation of commutators, which makes 
them suitable for the efficient propagation of systems with a large number of degrees 
of freedom. We present an optimized fourth order propagator and demonstrate its 
efficiency in comparison to the direct Rungc-Kutta computation. As an illustrative 
example we consider the parametrically driven dissipativc Dicke model, for which we 
calculate the periodic steady state and the optical emission spectrum. 
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1. Introduction 

The outcome of many experimental measurements is well described by linear response 
theory for situations close to thermal equilibrium. Other experiments, predominantly 
those dealing with small quantum systems in strong external fields, require a full 
non-equilibrium description. One example is cavity quantum electrodynamics pQ, and 
generally finite quantum systems in radiation fields. While the interaction of a single 
atom or an atomic ensemble with the quantized cavity field is weak, transitions between 
atomic levels can be induced with strong, classical laser fields. Through cavity losses and 
spontaneous emission the energy input from the external pumping dissipates. Atoms in 
a cavity are open quantum systems far from equilibrium. 

Such situations are described either by the Schrodinger equation 



with a time-dependent Hamilton operator H(t) if we neglect dissipation, or more 
generally by a master equation 



with a time-dependent Liouville operator C(t), e.g. one of Lindblad-type which includes 
dissipation in the Markovian approximation [2]. In addition to single-time expectation 
values, which provide the basic information from time propagation of the wave function 
or density matrix p(t), one is interested in many-time correlations functions that 
yield optical spectra or information about the coherence or statistical properties of the 
emitted light 

Since explicit solutions of linear differential equations with variable coefficients 
do not exist apart from simple situations, the above equations fall into the domain 
of numerical time-propagation. The topic of the present paper is the application 
of commutator- free propagators based on the Magnus series [5]. The Magnus series 
arises in the context of differential equations on Lie groups, where it allows, among 
many other things, for the systematic construction of high-order approximations to the 
propagator j6]. Commutator-free exponential time-propagators (CFETs) avoid the use 
of commutators that appear in the Magnus series. They provide an efficient and accurate 
algorithm for numerical time propagation [7j, which we discussed for the Schrodinger 
equation in reference [Sj. Here, we concentrate instead on master equations for open 
quantum systems. Although the present application lies outside of the principal Lie 
group setting, we feel that the numerical results presented here are promising enough to 
warrant closer inspection. The application to the parametrically driven dissipative Dicke 
model in section [5] gives an indication of the potential of this approach in non-trivial 
situations. 

The paper is organized as follows. In sections |2] and [3] we discuss the basic numerical 
problem and its principal solution through the Magnus series. The commutator-free 
exponential time-propagators are introduced as a more practical solution in section HJ 
and an optimized 4th-order propagator is given in 14.21 After a demonstration of their 
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usage with the example of a spin in a magnetic field in section I4T31 we turn to a discussion 
of the parametrically driven dissipative Dicke model in section |5l before we summarize 
in section El 

2. The numerical problem 

To describe the basic numerical problem we consider the Schrodinger equation ([1]). The 
standard approach obtains the wave function \ip(t)) through time stepping. The middle- 
point approximation 



allows for propagation of \ip(t)) over a short time interval [t, t+St]. Repeated application 
of equation ([3]) gives \ip(t + T)) starting from \ip{t)) with N = T/5t time steps. As 
detailed later, straightforward expansion of the exponential shows that the error of one 
time step with equation ([3]) is oc (St) 3 , such that the total error oc N(St) 3 = T(5t) 2 for 
propagation time T scales as (St) 2 . Conversely, the achieved accuracy scales as (St)~ 2 , 
which we can write symbolically as error oc effort -2 . 

As a second-order method the middle-point approximation is not efficient and 
requires small St even for low accuracy demands. If we ask for a better scheme we 
should note that the approximation has two independent sources of error. The 
genuine error in the situation of a time-dependent Hamilton operator arises from the 
replacement of H(t) by the constant H(t + 6t/2), and depends mainly on the rate 
of change of H(t). In addition, the numerical computation of an operator or matrix 
exponential exp[A] involves an error determined by the spread of eigenvalues of A. In 
equation ([3]) it is roughly proportionally to 5t and the norm \\H(t + 5t/2)\\. 

Often the rate of change of H(t), e.g. set by an external field frequency, is smaller 
than the largest eigenvalues of the Hamilton operator corresponding to highly excited 
states. Then the total error is dominated by the computation of the exponential in 
equation ([3]). Very small time steps St and correspondingly large effort are required 
even if H(t) changes slowly. 

This observation explains why the use of general algorithms for the solution of 
ordinary differential equations (ODEs), e.g. the standard 4th-order Runge-Kutta (RK4) 
procedure [9], cannot be recommended unreservedly for the Schrodinger equation. The 
problem of such (explicit) ODE solvers is that they provide only a poor approximation 
of the exponential in equation ([3]) and are inefficient already for constant H. 

For example, the RK4 procedure approximates the exponential by the 5 terms 
exp[A] « 1 + A + A 2 fl + A 3 /6 + A 4 /2A of the Taylor series of e x . The problem 
is that the Taylor series is not a good approximation unless \x\ is very small. This 
effect is shown in panel (a) in figure [U where it is compared to an approximation using 
Chebyshev polynomials (of the first kind) [10J. In this example, the five term Chebyshev 
approximation is 16 times more accurate than the five term Taylor series. For the 4th 



\vj>(t + St)) « exp - i St H{t + St/2) \i/>(t)) 



(3) 




Figure 1. Left panel (a): 4th-order Taylor (green curve) and Chebyshev 
approximation (red curve) of Re exp(i) = cos(i) on the interval [—2,2]. The dashed 
lines gives the error of both approximations. The maximal error is 7.8 x 10~ 2 for 
the Taylor approximation, which loses accuracy at the boundaries of the interval, 
versus 4.7 x 10~ 3 for the Chebyshev approximation. Right panel (b): Error e (see 
equation (0|) for the calculation of cxp[— \Ht] with the diagonal 10 x 10 matrix with 
eigenvalues H nn = n and t = ir/5. We compare the 4th-ordcr Runge-Kutta procedure 
(RK4) with the use of the 4th-order Chebyshev approximation in a time-stepping 
scheme (Cheb 4 ), and with a single propagation step using N terms of the Taylor series 
(Tay) or of the Chebyshev approximation (Cheb). In time-stepping the error decays 
as a power (here cx ./V -1 / 4 ) of the effort, while full computation of the exponential in 
a single step achieves much quicker error reduction. 



order approximation error oc effort 4 , this implies that the efficiency is increased by a 
factor I6 1 / 4 = 2. 

In panel (b) in figure Q] we compare the error-effort relation of the RK4 procedure 
to that of the Chebyshev approximation for the calculation of exp[— iStH], where H is 
a diagonal matrix with entries H nn = n as for a harmonic oscillator. Here, and also 
in later examples (figures |2] and |3]) , we give the error between an exact and numerical 
matrix A e / n as the maximal difference of matrix elements 

e = max \A%-A%\ . (4) 

The Chebyshev approximation is clearly superior already for a small St. It becomes 
even better with increasing eigenvalue spread or time-step \St\. 

This reasoning motivates the replacement of ODE solvers by techniques which take 
advantage of the linearity of the Schrodinger or master equations. The calculation of 
a matrix exponential is better accomplished with specialized algorithms, such as split- 
operator methods [TT], or the Krylov [12j[T3] or Chebyshev technique [H] in the case 
of large sparse matrices. They allow for efficient propagation with time-independent 
Hamilton operators. Equipped with such algorithms it remains to improve on the 
genuine error oc (St) 2 involved in equation fl3]) when turning to time-dependent Hamilton 
operators. 
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3. Magnus propagators 

The importance of accurate evaluation of exponentials for the Schrodinger equation is 
related to the fact that the exponential maps the hermitian Hamilton operator H onto 
the unitary propagator exp[— \tH\. Many differential equations involve a Lie algebra 
(here: of hermitian Hamiltonians) and a Lie group (here: of unitary propagators) in 
this way. The idea behind "geometric numerical integration" of ODEs [15] is that also 
an approximate propagator should stay in the respective Lie group. 
Let us consider general linear differential equations 

±{t) = A{t)x{t) , (5) 

where x(t) is a vector and the coefficient matrix A(t) is time-dependent. The formal 
solution of equation fl5]) is provided by the propagator U (t) that gives 

x (t) = U(t)x(0) (6) 

for all x(0). 

For a scalar equation x = a(t)x, the propagator is obtained through integration 
U (t) = exp[J Q * a(T)dr]. For operators or matrices [A(ti), A(t 2 )] 7^ is possible, such that 
this expression does not generalize. We can still read this expression as an approximation 

U{t) « exp [ [ A(t )dr 1 . (7) 







For small t = 5t, this is nothing else than the approximation fl3]), if the integral over r 
is approximated (also with error {St) 3 ) using the middle-point value A(5t/2). Although 
equation ([Tj) involves a finite, maybe large, error its exponential form guarantees that 
the approximate propagator lies in the Lie group. The question is whether we can 
improve on the (St) 3 scaling of the error and preserve the exponential form. 
An affirmative answer is given by the Magnus series [SJE], which gives 

U(t) = exp [ f A( T )dr + Q 2 (t) + Q 3 {t) + . . . ] (8) 
L Jo J 

as an exponential of Lie algebra elements £l n (t) and provides a systematic scheme for 

their construction. 

The first term in (jSj) is the term known from the scalar case. The non-commutativity 
of A(t) is accounted for by correction terms Q n (t), for n > 2. The term Q n (t) is given 
by a time-ordered integral of n-fold nested commutators of A{ji) and can be obtained 
through a recursive calculation. The first two terms are 

Mt) = I I dn r dr 2 [A(n), A(r 2 )} (9) 
^ Jo Jo 

and 

n 3 (t) = ^^dr 1 £dr 2 £dT 3 [A(T 1 ) J [A(T 2 ),A(T 3 )]} + [[A(n), A(r 2 )}, A(r 3 )\ . (10) 
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Explicit expressions for higher-order terms become quickly unwieldy. Importantly, by 
building terms from commutators of Lie algebra elements Afc), every Q n (t) stays in 
the Lie algebra. 

The Magnus series solves two problems. On the one hand, it preserves the 
Lie group structure of a differential equation. For the Schrodinger equation, where 
A(t) = —iH(t), the propagator U(t) is unitary as the exponential of an anti-hermitian 
matrix. Furthermore, unitarity of U(t) is preserved for any truncation of the Magnus 
series. On the other hand, since the term Q n (t) involves an n-fold integration over time, 
its size scales as {St) n . Working with a truncated series including terms Vt n for n < N 
only, the error of the obtained approximate propagator itself scales as (5t) N+1 . We can 
thus improve systematically on the middle-point approximation (j3J) by including more 
terms from the Magnus series. 

Unfortunately, the Magnus series does not solve the practical problem of 
finding an efficient numerical time-propagation algorithm. Computation of the 
nested commutators and multiple integrals is difficult to implement and consumes 
computational resources. Fortunately, there is a simpler and more convenient way. 



4. Commutator-free exponential time-propagators 

The use of commutator-free exponential time-propagators (CFETs) has been discussed 
in references [HIEIEEE]. They are, basically, a reformulation of the Magnus series that 
avoids integrals and commutators and gives the propagator as a product of exponentials 
of simple linear combinations of A(t). 

The simplest CFET is the middle-point approximation itself. A 4th-order 
CFET, where the error scales as (5t) 4 , was introduced in [7UT6] . It gives the approximate 
propagator as the product of two exponentials 

U CFET (St) = exp [st ( 9l AW + g 2 A^) ] exp [st (g 2 A^ + 9l A®) ] , (11) 

which involve a linear combination of A(t) specified by the coefficients 

3-2V3 3 + 2V3 

91 = ^2— ' 92 = ^2~ ' (12) 
and uses only the values 

A« = A[ Xl 5t) , A^ = A[x 2 8t) (13) 

of A(t) evaluated at two points in [0, St] given by 

1 \/3 1 V3 

xi = , x 2 = -H • (14) 

2 6 2 6 V ; 

This expression is the simplest non-trivial CFET. A better, optimized, 4th-order 
CFET is presented below in equation ( 1221) . Higher-order CFETs can be constructed, 
and the freedom in the choice of coefficients can be exploited for their optimization, i.e. 
the minimization of the error. The construction of CFETs is rooted in the theory of 
abstract free Lie algebras that underlies the Magnus series. Its description is beyond 



Numerical time propagation of quantum systems in radiation fields 



7 



the scope of this paper, and we refer the reader to reference [7] and our reference [5] for 
details. Here, we proceed in the opposite way and give a direct check of the validity of 
equation f TTTT) that avoids most of the language of free Lie algebras. 



4-1. Direct validation of the 4th- order CFET 

The principle idea is to combine the two exponentials in equation (ITT]) with the Baker- 
Campbell-Hausdorff (BCH) formula exp[X] exp[F] = exp[X + Y + [X, Y]/2 + ...] and 
compare the resulting expression with the original Magnus series. Let us begin with the 
Taylor series 

A(t) =A, + A 2 t + A 3 t 2 + A 4 t 3 + 0(t 4 ) (15) 

of A(t), in the vicinity of t = 0. For the 4th-order CFET, only terms Ay, . . . , A A have 
to be considered. 

We insert the Taylor series in the Magnus series (jHJ), and keep the first 
three terms to Q A (t). The terms Q n (t) for n > 5 give contributions of order 
(St) 5 and higher. A simple counting of indices shows that only the seven terms 
Ai, Aii A3, A4, [Ay, A 2 \, [Ai, A 3 ], [Ai, [Ai, A-2]] can contribute in fourth order. The 
Magnus series thus gives the propagator 



U(St) = exp 



6tA l + { -^A 2 + (6t) 3 ^A 3 - ±[A X ,A 2 
+ (tty(\A 4 -±[A u A 3 ])+0((6t) 5 



(16) 



Note that the commutator [Ax, [Al,^]] does not contribute. This expression is the 
exact reference for comparison. 

It is easy to see that the middle-point approximation ([3]) is correct to second order 
(St) 2 : The terms StAi + ((St) 2 /2)A 2 are reproduced exactly, but the commutator [Ax, A 2 ] 
in the third order term is missing. 

For the 4th order CFET from equation (fTTl) . it is 

A^ = Ax + 5tx k A 2 + (5tx k ) 2 A 3 + (6tx k ) 3 A 4 + 0((5t) 4 ) (17) 
for k = 1,2, which inserted gives 

f/ CFET (5t) = exp 6t( gi + g 2 )Ax + (5t) 2 (g lXl + ^2)^2 (18) 

+ (5tf( gi x\ + g 2 x\)A 3 + (bt)\g x x\ + g 2 x\)A A + 0((5tf) 

x exp St(g 2 + gx)Ax + (St) 2 (g 2 xx + #1X2)^2 

+ (5tf(g 2 x\ + gxx 2 2 )A 3 + (5t)\g 2 x\ + gxx\)A A + 0((5tf) 
We now use the BCH formula 

exp \x + Y + l[X, Y] + hx, [X, Y}} - l[y, [X, Y}} + . . .1 (19) 



e x e Y 
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to combine the two exponentials. Only the three commutators shown in equation (1191) 
need to be evaluated, the following commutators in the BCH formula contribute terms 
of order (5t) 5 or higher. We then obtain the expression 

UcFET(St) = exp [sfaAi + {5t) 2 i 2 A 2 + (5t) 3 (^ 3 A 3 + A 2 ]) 

+ {5tY{i A A A + X2 [A l) A,} + x,[A l) [A l) A 2 ]]) +0((5t) 5 )} , (20) 

which allows for direct comparison with the Magnus series in equation ( 1T6l) . We 
immediately recognize the seven terms and commutators from there. 

Comparison of the coefficients which we evaluated with the BCH formula, to 
the prefactors in ( 1T6l) gives the conditions 



£l = 


2 9l + 2g 2 = l, 






(21a) 


£2 = 


(9i +92){xi +x 2 ) = - , 






(216) 


£3 = 


(9i + 92)(xl + x 2 2 ) = - , 






(21c) 


£4 = 


(gi + 92){4 + 4) = \ » 






(21 d) 


Xi = 


^(^1 +^2) ((5-2^1 +9ix 2 ) - 


(gixi + 92x 2 )^ 


1 

~ ~12 ' 


(21 e) 


X2 = 


^(9i +92) {[92x\ + 9ix 2 2 ) - 


{gix\ + g 2 xl)^ 


1 

~ "12 ' 


(21/) 


X3 = 


Y^(9i + 92) 2 (^{92X1 + gix 2 ) 


- (92X1 + gix 2 ) 


)-o. 


(21^) 



To complete our check we insert x%, x 2 from equation (|T4l) and ^1, ^2 from (fl2"|) and find 
that all the seven conditions are satisfied. Note that condition ( |21gP , and also ( 12 Id) 



and ( |21J , are redundant. 

For the construction of a CFET, this process has to be reversed. We start from 
an ansatz with a product of exponentials, derive the relevant conditions for a CFET 
of required order, and solve the resulting polynomials equations for possible coefficient 
values. Several adjustments of the direct calculation done here simplify bookkeeping, 
and reveal underlying structures which reduce the number of conditions. Still, the 
construction of higher-order CFETs is involved and not entirely free of brute force 
computations. We refer the reader to reference [8] to get an impression, where CFETs 
up to order 8 are presented. The usage of a given CFET, however, is plain and simple. 

4-2. Optimized 4th order CFET 

For the application to dissipative systems we recommend here the use of an optimized 
4th-order CFET, which is equation (43) in our reference [8]. Extending the simpler 
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expression (ITTj) it gives the approximate propagator as the product of three exponentials 
u cFwr( 5t ) = ex P 



x exp 
x exp 



St (M W + 92A^ + g 3 A®) 
St (g 4 AW + 95 A^ + g 4 A {3) ) 
St {g 3 AU + g^*) + giA ^) 



(22) 



where 



with 



and 



A^ = A[x x 5t] , A® = A[x 2 St] , A® = A[x 3 St] 



fji 



94 



1 

2 ~ 

37 
240 




%2 



1 

2 



10 

87 



fj2 



X 3 
1 

30 



1 

2 + 

</3 




(23) 



(24) 



37 10 
240 + 87 



11 



360 



23 

95 = 45 



(25) 



The error of this CFET scales again as {St) 5 , but the prefactor in front of the error 
term is considerably smaller than for the CFET (fi~T|) . The reduction outweighs the 
increase of effort using three instead of two exponentials. 

In contrast to the original Magnus series usage of the CFET (1221) is compellingly 
easy. Only linear combinations of A(t) evaluated at three different points in the interval 
[0, St] need to be formed. All commutators and integrations have been removed from 
the expression. Of course, we assume the existence of an algorithm for the computation 
of matrix exponentials. 

Let us stress the advantage of easy usage with the even simpler formulation that is 
obtained for an A(t) = B + f(t)C (it is easily generalized to include more terms). In 
this case equation (122]) can be written as 



^cfet(^) = exp 
with the time steps 

Stx = 



Sh (B + fiC) 



exp 



St 2 (B + f 2 C) 



exp 



5h (B + f 3 C) 



11 P 
— St 
40 



Su = —St 



20 



and the coefficients fi, f 2 , from 



( fA 




( h h 2 h 3 \ 




( fixiSt) \ 






/I5 /I4 




f(x 2 St) 


\ h ) 




\ h h 2 hi J 




\ f{x 3 St) j 



(26) 
(27) 

(28) 



using 



hi 
h 4 



37 400 

66 ~ 957 
11 

~ 162 ' 



ho = - 



33 



37 
66 



400 
957 



92 
81 



(29) 
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Through the CFET the full propagation with a time-dependent term f(t)C is 
replaced (approximately) by piecewise propagation with a constant term /jC. The 
choice of the coefficients fi, fz, fz according to equation ( )28|) guarantees that the error 
of this approximation scales as (St) 5 . This is the signature of geometric integration [T5] : 
Figuratively speaking, instead of moving along a curve in the Lie group we move 
repeatedly along short straight lines, the direction of which is given by the linear 
combinations in equation ( 122]) or (126]) . 

Note that the time steps St±, St2 are positive such that propagation proceeds in the 
forward direction. This is important for dissipative systems where a negative 5ti would 
push eigenvalues of C(t) into the right half complex plane, which leads to exponentially 
growing terms and corresponding numerical instabilities. For this reason we restrict 
ourselves to 4th-order CFETs here. 

4-3. Exemplary application of the optimized J^th-order CFET 

Let us apply the 4th-order CFET ( )22l) to a simple example and compare with the RK4 
procedure which, as we claimed, should be less efficient because it does not properly 
compute exponentials. We have to stress that the efficiency of CFETs depends on a good 
algorithm for the computation of matrix exponentials. Otherwise, when the additional 
computational overhead involved exceeds the savings achieved with a large time-step 5t, 
the simple Runge-Kutta procedure is more efficient. 

Good algorithm for the symmetric case (A = ±A, e.g. for a hermitian 
Hamiltonian) are the Chebyshev, Krylov and split-operator techniques mentioned 
before. For the unsymmetric case ([A, A'] ^ 0) encountered for dissipative systems, all 
techniques meet problems which are only partially solved. After suitable modifications 
of the standard procedure the Chebyshev technique behaves most favourably. The 
exploration of this point has to be left for a future publication, here we take the virtues 
of the Chebyshev technique for granted. 

4-3.1. First example: Spin in a magnetic field As an example for a non-dissipative 
system consider a spin (length j) in a rotating magnetic field, with Hamilton operator 

H(t) = 2AJ Z + 2V cos 2utJ x + 2V sin 2utJ y . (30) 

The time-evolution of the wave function can be determined exactly after transformation 
with exp[iutJ z ] to the rotating frame (see below). 

In figure [2] we plot the error-effort relation for the optimized CFET ( 122]) and the 
RK4 procedure for a typical set of parameters (the behavior for other parameters is 
identical). The error e is determined as in equation (111). For the effort Nh we count 
the number of evaluations of matrix-vector multiplications with the Hamilton operator, 
which is generally the most time-consuming step. The matrix exponentials needed for 
the CFET are calculated with the Chebyshev technique to machine precision. 

We see that the RK4 procedure requires lesser effort for low accuracy only, but 
use of the CFET becomes quickly advantageous as the spin length or propagation time 
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Figure 2. Error e between exact and numerical density matrices p e / n (t) (see 
equation (QJ) versus the effort Njj (number of matrix- vector multiplications with H) 
for A = V = OJ = 1.0 and j = 1/2 (left panel (a)) or j = 10 (right panel (b)). The black 
curve gives the results using CFETs with a Chebyshev evaluation of the exponential, 
the red curve using the RK4 procedure. Curves are shown for propagation time t = 10 
and t = 100. The grey dashed lines indicate the reduction of the RK4 effort by a factor 
1/2 or 1/4 achieved by the CFETs. 



increases. For j = 10 and t = 100 in panel (b), the CFET is more efficient for error 
goals less than 1%, with an a efficiency gain of a factor 2 ... 4. 

4-3.2. Second example: Driven dissipative two-level system We keep the spin in the 
rotating magnetic field as an example and include dissipation. With dissipation, its time 
evolution is described by a master equation ((2]) for the spin density matrix p(t). 

In the Lindblad formalism, the Liouville operator £(t) = Cnif) + C D is the sum of 
two or more terms with different meaning [2]. The first term 

C H (t)p = -i[H(t),p] (31) 

contains the Hamilton operator H(t) and will be time-dependent. The second and 
further terms have the form 

V[A\p = lApti - A ] A P - P A ] A . (32) 

They introduce eigenvalues of C(t) with finite (negative) real part and thus describe 
dissipation. Within the Lindblad formalism the form of T>[A] guarantees that the 
structural properties of the density matrix - hermiticity, normalization, positive semi- 
definiteness - are strictly preserved. Note that the numerical time propagation scheme 
does not depend on the precise form of C(t), as long as the master equation remains 
linear and local in time. 

For the driven spin, dissipation is included through the Lindblad term 

V[Jjp = 2J„pJ + - J + J^p - pJ + J_ , (33) 

and the full Liouville operator is 

C = C[H] +tZ>[J_] (34) 
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Figure 3. Left panel (a): (Jz(t)) of the driven dissipative spin j = 1/2 from 
equation with A = V = oj = 1.0 and 7 = 0.01. Shown is the envelope function, 
suppressing the fast spin oscillations with frequency A. The magnetic field frequency 
u)(t) grows linearly from to 2 during propagation (green dashed curve). The red 
curve gives the steady state result from equation (|36[) to the respective frequency ui(i). 
Right panel (b): Similar to figure [2j the error e vs. effort Njj for A = V = u> = 1.0, 
j = 1/2 and finite dissipation 7 = 0.01. 



with the dissipation rate 7 > 0. 

For j = 1/2, the exact solution of this problem is possible with a transformation 
p{t) = exp[io;t J z ]p(t) exp[— iut J z ] to the rotating frame, which gives a time-independent 
Hamilton operator H = 2(A — oj) J z + 2VJ X . Note that the transformation leaves J z 
invariant. 

The stationary state in the rotating frame, corresponding to the eigenvalue zero of 
the transformed Liouville operator C for 7 > 0, is 

_ 1 / V 2 -(2(A-w)+i 7 )7 \ 

Po ° 4(A - ujf + 7 2 + 2V 2 ^ _( 2 (A - U ) - \rj)V 4(A - uf + 7 2 + V* J ' 
In particular, (J z {t)) converges for t — > 00, with constant value 

V 2 1 

{Jz)co = 4(A- W ) 2 + 7 2 + 2V 2 " 2 ' (36) 

In panel (a) in figure [3] we plot (J z (t)) starting from the initial state with 
(Jj(0)} = +1/2. Transient oscillations decay as a result of finite dissipation 7 > 0. The 
value of (J z (t)) in the quasi-equilibrium state depends on the field frequency u, which 
we increase slowly during time-propagation. (J z {t)) follows closely the value ( J z )oo from 
equation (1361) . with a short delay, and we identify the resonance at uit) — A. 

At resonance u = A, we get simple expressions for the remaining three eigenvectors 
of C. The non-zero eigenvalues are Ai = —7, A2 = —(37 + 0/2, an d A3 = —(37 — 0/2, 
with corresponding eigenvectors 



Pi = h n h (37) 
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1 _i( 7 _£)/(4V) 

P2 ~ { i(7-0/(4V) -1 I' (3 ' 

1 -i(7 + 0/(4V) 

P3 -' Ki+o/m -i (39) 



We have introduced the abbreviation £ = \Ay 2 — 16V 2 . 

Starting from the initial state with (J^(0)) = |, we obtain the time evolution of 
from the decomposition 

/ 1 \ V 2 + 7 2 7 3 + 5>yV 2 

P(0) = ^ Q Q j=P^+ 2{2V2 - (P2 + P3) + 2e(2y2 + 72) (P2 - P3) (40) 

of the density matrix. In the underdamped case 7 < AV, we have 

-7 2 

(Ut)) - 



2( 7 2 + 2T/ 2 ) 



V 2 + 7 2 ) coscDt - ( 7 3 + 5 7 ^ 2 ) e" (3/2ht , 



2V 2 + 7 2 V 2o) 



where we write a) = (l/2)\JlQV 2 — 7 2 for the frequency of the transient contribution. 
These expressions allow for comparison with numerical results. 

For the numerical solution of this problem, we do not transform the problem but 
keep the time- dependence of H(t) explicitly. In panel (b) in figure [3] we plot the error- 
effort relation as in figure [2j The relation between the CFET and the RK4 procedure 
is similar to the non-dissipative case, and we recognize the factor 1/2 of error reduction 
for j = 1/2. The advantage of the CFET is however not quite as distinct as for the 
dissipation-free case in figure El panel (a). 



5. The parametrically driven dissipative Dicke model 

The previous examples serve as a benchmark for the CFET approach. We can now 
demonstrate its usefulness for a less academic case, the parametrically driven dissipative 
Dicke model. Optical properties of the Rabi case j = 1/2, corresponding to a single 
qubit, have been explored in reference [17] , and quantum phase transitions in the 
parametrically driven Dicke model without dissipation are studied in [18]. 
The Hamilton operator of the Dicke model |19j . 

H(t) = &J z + VLa)a + \(t)(a + a))J x , (42) 

describes an ensemble of two-level atoms (transition energy A) as a pseudo-spin of 
length j, which couples to the cavity field (frequency Q). We assume a time-dependent 
interaction constant 



\(t) = Aq + 5 \ cos Upt , 



(43) 
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Figure 4. Comparison of the recommended (equation (|22|l) and simpler 
(equation (fTTj) ) CFET with the RK4 procedure and the middle-point approximation © 
in application to the parametrically driven dissipative Dicke model. The Dicke model, 
with parameters A = fi = 1, n = 0.01, An = 1, 6X = 0.5, u p = 2, is propagated 
over 30 modulation periods, i.c for < t < 30(2tt/u! p ). As in figure [2] we show the 
error e versus the effort Nh for all propagation schemes. The error is determined from 
comparison of the numerical density matrix and the reference solution obtained in the 
limit Nh — > oo. The CFETs are used in combination with a Chebyshev evaluation 
of the matrix exponential. Left panel: For j = 1/2, the CFET (f22|) is 4 times more 
efficient than the RK4 procedure. Right panel: For j = 5, the CFET (|22|) is 8 times 
more efficient than the RK4 procedure. 



and consider dissipation only through cavity losses described by the term T>[a] (see 
equation (152]) ) but neglect spontaneous emission (e.g. described by D[J_]). The 
Liouville operator is 

C = £ H (t) + nV[a] , (44) 

with loss rate k > 0. In contrast to the standard quantum optical treatment in rotating 
wave approximation, it is not possible to eliminate the explicit time-dependence of H(t) 
through a transformation to the rotating frame. 

Because the rotating wave approximation is not applicable, physical properties 
of the parametrically driven dissipative Dicke model have to be extracted from time 
propagation of the density matrix of the joint atom-photon system. The number of 
entries of the density matrix, which grows as oc (2j + l) 2 , becomes large already for 
moderate pseudo-spin length j. In addition, highly excited states contribute to the 
dynamics if j grows. Therefore, as we discussed in section [21 the advantage of CFETs 
over general ODE solvers such as the RK4 procedure will be pronounced for this more 
complex example. 

In figure H] we compare the recommended CFET from equation (l2"2l) to the RK4 
procedure and the naive middle-point approximation ([3]). We see that already for 
j = 5 we can easily reduce the numerical effort by a factor eight if we use CFETs. 
The reduction is achieved independently of the intended error e. The middle-point 
approximation, which is only a second order scheme, is not able to compete even with 
the RK4 procedure. This clearly supports our recommendation for the CFET (122j) . and 
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Figure 5. Panel (a): Time-dependence of I z {t) in the parametrically driven dissipative 
Dicke model for j = 5, with a slow linear change of the modulation frequency tu p 
from 1.5 to 2.5. The parameters are = A = 1, « = 5 x 10~ 3 , Xq = 0.02, and 
SX/Xo = 2.5 x 10~ 2 . Shown is the envelope function of I z (t) instead of the fast 
oscillations with frequency A. Panel (b): Lowest energy levels of the Dicke model 
at weak coupling and = A. In the Rabi case j = 1/2 only two doubly excited 
states with energies £/2,± exist, for j > 1/2 an additional third state with energy £2,0 
appears. Transitions from the ground state change the number of excitations by two, 
and follow the vertical arrows. The matrix element for the transition with energy E^o 
vanishes in lowest order perturbation theory. 

extends the positive results from [5] to dissipative systems. Note that the CFET (I22p is 
better (by 50%) than the simpler CFET (11 ip although it requires computation of three 
instead of two matrix exponentials. 

Note again that the advantage of CFETs stems from the fact that the propagator 
is approximated as an product of exponentials of the Hamilton or Liouville operator. 
We know that this strategy is favourable for non-dissipative systems because it respects 
the unitary geometry of the Schrodinger equation [6l[T5], but apparently it works well 
also for density matrix propagation in driven dissipative systems. Conversely, CFETs 
rely on good computation of matrix exponentials. As mentioned at the beginning of 
section 14.31 we use the Chebyshev technique here because it proved to be more efficient 
and reliable than a Krylov computation for non-hermitian matrices. 

5.1. Steady state resonances 

In panel (a) in figure |5] we show the time evolution of the initial state p(0) = 
where = \—j/2) (g) |vac) is the state with no atomic or field excitations. We plot the 
population inversion 

h{t) = \ + -U z {t)) (45) 

for a linear variation of the modulation frequency u p from 1.5 to 2.5. Transient 
oscillations are observed for t < 500, before two resonances evolve at a time 
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corresponding to w p ~ 2 ± 0.12. Beyond the resonances, I z (t) decays again to a small 
value with weak oscillations. 

The energy level diagram in panel (b) in figure [5] explains the appearance of 
resonances in panel (a) through transitions between the states of the Jaynes-Cummings 
ladder [20]. 

The eigenstates of the zero coupling Hamiltonian are the J z , a^a eigenstates \m, n), 
with m + j atomic excitations and n photons. For Q = A the states \m + k, n — k), for 
several integer k, are degenerate with energy (m + n)Q. At weak coupling Ao <C ^, A, 
degenerate states are split oc A by the atom-field coupling. Counting energies relative 
to the energy of the lowest state |— j, 0), the energies of the two singly excited states 
(|-j + 1, 0) ± \-j, 1))/V2 are given by E lt ± = Q ± VW^o- 

A weak modulation of X(t) introduces transitions between | —j, 0) and the doubly 
excited states (vertical arrows in panel (b)). Note that the splitting of energy levels 
in the diagram is determined by the co- rotating terms J_a\ J + a in the coupling term 
J x (a + a^), which preserve the number of excitations in the sense of the rotating wave 
approximation, but the transitions arise from the counter-rotating terms J+a', J_a and 
change the number of excitations by two. In the Rabi case j = 1/2, the two doubly 
excited states (| 1/2, 1) ± | — 1/2, 2>)/>/2 have energy E 2 ± = 2f2 ± y/2X . Resonances are 
expected at these energies [T7] . 

In the Dicke case j > 1/2, the splitting of the three degenerate states \—j + 2,0), 
\—j + 1, 1), 2) has to be calculated. This gives the energies E 2j ± = 2Q ± s/Ej — 2A 
for the odd/even parity combination, reproducing the j = 1/2 result. The third state 
with unshifted energy i^.o = 2fi is a linear combination of \—j + 2,0), |— j, 2) only. It 
does not couple to 0) through the counter-rotating terms, and the transition to this 
state is forbidden in leading order of perturbation theory. Therefore, we also expect 
only two resonances in the Dicke case. For the parameters from figure [5l they occur at 
E 2 ,± = 2 ± 0.02 V38 ^2± 0.1233. 

To identify the resonances numerically we propagate the system with fixed u p until 
the periodic steady state is reached. Then we calculate the quantities 

h= / L(t')dt\ N b = / N b {t')dt' (46) 

Jt Jt 

averaged over one modulation period 2n/uj p . Here, 

N b (t) = {a)(t)a(t)) . (47) 

is the number of cavity bosons. 

In figure [6] the quantities I z , N b are shown as a function of ui p . We recognize the 
two resonances u p ~ E 2 ,±, which are broadened due to the cavity losses oc k. For 
the calculation we used the optimized 4th-order CFET (|22|) together with a Chebyshev 
computation of the exponential. 
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Figure 6. Population inversion I z (black curve), number of cavity bosons JV& (red 
curve), and total emission S t ot (dashed red curve) for the driven dissipative Dicke 
model, as a function of the modulation frequency ui p with parameters il = A = 1, 
K = 5 x 10~ 3 , A = 0.02, <5A/A = 2.5 x 10~ 2 as in figure [51 The quantities arc 
averaged over one modulation period 2tt/ui p . Panel (a): Results for the Rabi case 
j = 1/2. Panel (b): Results for j = 5. The grey dashed lines indicates the resonances 



at LOr 



E- 



2,±- 



5.2. Emission spectrum 

To study the optical properties of this system we compute the cavity emission spectrum 
S(u). It is obtained as the Fourier transform 

1 r°° 

S(u) = - Re / S(r)e- luJt dr (48) 
of the correlation function 

pt-\-2n /uj p 

S(r) = / (a\t + T)a{t'))dt' , (49) 



which we calculate with the quantum regression theorem [5] through time propagation 
of the operator ap(t) (for r > 0). The correlation function involves the average over one 
modulation period [t,t + 2tx /u p ] for large t, i.e. in the periodic steady state. 
We include in figure E] the total emission 



S tot = / S(u)du>, (50) 



o 



which is given by the integral over positive u in accordance with the fact that emission 
of a (real) photon can only decrease the energy. We note the normalization 

S(u)du = S(t = 0) = N b . (51) 

> 

We see that S tot ~ Nf, close to resonance, when emission is strong. Away from resonance 
Stot drops below iVf,, since iV& counts also bound photons that do not contribute to 
emission. However, Stot remains finite as a consequence of the Markovian approximation 
used here for the dissipative term [17] . 

The emission spectrum for the Rabi case j = 1/2 is shown in figures [71 [8] and for the 
Dicke case with j = 5 in figures [91 [TUl For weak coupling, i.e. for |A(t)| {A,f2}, the 
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Figure 7. Emission spectrum S(ui) for the Rabi case j = 1/2, at (panel (b)) and 
close to (panels (c), (d)) the higher resonance uj p = £2,+- It is oj p = £»2,+ — 0.0080 in 
panel (c) and lo p = E^,+ + 0.004O in panel (d). The remaining parameters are as in 
figure [6] The energy level diagram in panel (a) follows figure The four transitions 
are marked by vertical arrows and indicated by corresponding numbers in the other 
panels. The transition energies are given in equation (|52j) . The spectrum from panel 
(b) is included in panels (c), (d) as the grey filled curve. 
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Figure 8. Emission spectrum S(u>) for the Rabi case j = 1/2, at (panel (b)) and 
close to (panels (c), (d)) the lower resonance uj p = E-2.-- It is uj p = E2- — 0.0040 in 
panel (c) and uj. p = E 2 ,- + 0.008f2 in panel (d). The remaining parameters are as in 
figure El and the notation follows the previous figure [71 
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Figure 9. Emission spectrum S(u>) for j = 5, at (panel (b)) and close to (panels 
(c), (d)) the higher resonance w p = £2,+ - It is ui p = £2,+ — 0.04f2 in panel (c) and 
uj p = £2,+ + 0.04f2 in panel (d). The remaining parameters are as in figure [6j and 
the notation follows the previous figures The transition energies are given in 
equation ([53]) . The weak signals of the parity- forbidden transitions 4, 4' are displayed 
with magnification factor 20 in all panels. 



interpretation of the emission spectrum is again possible using the energy level diagram 
from figure [5j 

For Up ~ £2,+ in figure [TJ the higher of the two doubly excited states is populated. 
Since the operator a changes the number of excitations by one, an atom in this state 
can decay to the lowest state only through the intermediate singly excited states. Four 
transitions corresponding to the four vertical arrows in panel (a) can be identified in 
this situation. For the Rabi case j = 1/2, they are in order of increasing energy 

Wl = E x _ =Q-X ^ 0.98 , u 2 = E 2t+ - E h+ = Q + A (v / 2 - 1) « 1.01 , 

w 3 = .Ei )+ = fi + Ao«1.02, u A = E 2t+ -Ek- = tt + A (v / 2+ 1) ~ 1.05 . (52) 

The numerical values correspond to the parameters from figure and the transitions 
are marked correspondingly in panels (a), (b). 

In transitions 2 and 4 the doubly excited state decays. These transitions are 
shifted by u p — E 2j + away from resonance (transitions 2', 4' in panels (c) and (d)). 
In transitions 1 and 3 the singly excited states decay, and their energy does not depend 
on the modulation frequency. Note that transition 4 (dashed arrow) is parity forbidden 
in lowest order perturbation theory and gives only a weak signal. 

The opposite case u p ~ E 2 _ in figure [8]has an analogous interpretation that follows 
from the energy diagram in panel (a). Now the lower of the doubly excited states is 
populated, and the order of transitions is reversed, with transition 1 as the parity 
forbidden transition. 

For the Dicke case j > 1/2 in figures 191 fTUl we expect the same qualitative behavior 
as for j = 1/2 since the additional state with energy E 2 G does not participate in the 
transitions. For the higher resonance u p ~ E 2>+ in figure M, the situation corresponding 
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Figure 10. Emission spectrum S(ui) for j = 5, at (panel (b)) and close to (panels 



(c), (d)) the lower resonance uj p = E-i—- It is uj p 



E 2 



0.04J7 in panel (c) and 



uj p = £2,- + 0.04O in panel (d). The remaining parameters are as in figure [6j and 
the notation follows the previous figures [7H9] The weak signals of the parity- forbidden 
transitions 1, 1' are displayed with magnification factor 20 in all panels. 
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Figure 11. Peaks in the emission spectrum S(io) at stronger coupling, for the Rabi 
case j = 1/2 at the higher resonance lu p = i?2,+- Left panel: Energies w pC ak of the 
different peaks as a function of coupling Ao. The solid curves marked 1-4 correspond 
to the transitions in figure [7] at Ao = 0.02, the dashed curves show the position of 
additional peaks emerging at stronger coupling. Because the peaks in S(uj) have a 
finite width, curves in the panel can begin and end in an isolated point. Right panel: 
Peak height S^cjpoak) of the five highlighted curves in the left panel. The height of the 
parity forbidden transition 4 (cf. Panel (a) in figure [7]) is close to zero for all Aq. 



to figure the transitions in order of increasing energy are 

u 1 = Q- « 0.94 , u 2 = ft + A ( v / 8j - 2 - y/2j) w 1.06 , 

u 3 = n + \ oy /2j « 1.06 , w 4 = A + A (v / 8j - 2 + y/2j) « 1.19 . (53) 

The numerical values correspond to the parameters from figure 

As a difference to the Rabi case j = 1/2 we note that w 2 ~ w 3 , such that the two 
transitions cannot be distinguished in panel (b) because of the finite linewidth acquired 
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through cavity losses. If we change u p transition 2 is shifted but transition 3 remains 
fixed. This allows for the separate identification of both peaks in panels (c), (d). In the 
same way we can understand the opposite case u p ~ E 2 - in figure HOI which is similar 
to figure EJ 

For stronger coupling, additional peaks appear in the emission spectrum. The peak 
position and height is shown in Fig. [11] for the situation corresponding to figured Note 
that we adjust the modulation frequency u p = E 2t + for different Ao to remain close to 
the higher resonance. For Ao > 0.1 a fifth transition peak "(5)" becomes visible, and 
the previous interpretation of S(u) via the Jaynes-Cummings ladder breaks down. At 
least such situations require numerical time propagation because a simple perturbative 
interpretation is no longer possible. 

6. Summary and Outlook 

Numerical time propagation allows for the theoretical description of experimentally 
relevant non-equilibrium situations beyond the linear response regime. Such situations 
arise in particular if small systems such as atoms or molecules are manipulated by 
strong radiation fields. Important directions of research include the optical properties of 
ensembles showing collective behavior, such as polariton or exciton condensates j2TJ[22] . 
A characteristic optical signature, as of the spatial shape and energy distribution of the 
optical emission j23j|2l] or the coherence properties of the emitted light, can provide 
the proof of existence for a condensate. A fundamental theory of optical properties 
of collective phases of (quasi-) particles with finite lifetime requires the description of 
the driven open quantum system that is realized, e.g., by the excitonic condensate in a 
semiconductor [25] . Different but on the fundamental technical level related questions 
arise in the field of non-equilibrium transport problems [26] . 

Increasing complexity of the physical situation under study coincides with an 
increase of the computational effort, which underlines the need for powerful numerical 
algorithms. We argued here in favor of commutator-free exponential time-propagators 
as a convenient alternative to the original Magnus series. Shifting the focus of previous 
studies, where CFETs were shown to be well suited for the time propagation of driven 
systems without dissipation, we here applied CFETs to open quantum systems. Using 
the parametrically driven Dicke model as an illustrative example we calculated the 
optical emission spectrum with this technique. 

Conceptually, CFETs are recipes for the reduction of the original problem - 
the solution of the Schrodinger or master equations with a time-dependent Hamilton 
operator — to the computation of matrix exponentials. Because they replace the naive 
approximation of the time propagator by a single exponential per time step, as it is 
encoded in the second order middle-point approximation fl3]), with a more sophisticated 
combination of exponentials, higher-order CFETs significantly reduce the numerical 
effort. The particular appeal of CFETs is that they can be combined with any technique 
for the computation of matrix exponentials, for example the powerful Krylov (Arnoldi) 



Numerical time propagation of quantum systems in radiation fields 



22 



or Chebyshev techniques in the context of large sparse matrix computations. CFETs 
do not compete with such techniques, but instead serve the complementary purpose of 
achieving a favorable error-effort scaling also for equations with a time-dependent H(t). 
Therefore, CFETs should be of interest to anyone presently using Krylov (Arnoldi) or 
Chebyshev techniques in studies of driven quantum systems: It is straightforward and 
simple enough to add the CFET computational scheme from equations ( 122]) . (126]) to an 
existing program or implementation [27] [28]. 

In addition to the principal research directions mentioned above many open 
problems arise within the more restricted context of the present work. A notorious 
problem is the efficient evaluation of the matrix exponential for non-symmetric large 
sparse matrices, which is essentially a problem of polynomial approximation in the 
complex plane without precise knowledge of the approximation domain. One may also 
question the principal usage of CFETs for dissipative systems, where dynamical semi- 
groups replace the Lie group setting. Modifications of CFETs can be tailored to this 
situation and circumvent the negative time-step problem that occurs for methods beyond 
the 4th-order CFETs to which we restricted our present considerations. 

A more physical question concerns the use of the Lindblad formalism in the 
description of optical emission. The Markovian approximation is not entirely 
satisfactory here, since it cannot distinguish between energy increasing ( "virtual" ) and 
energy decreasing ("real") transitions in the emission process. For this reason, the 
total emission rate remains finite (of the order of the cavity loss rate k) even without 
external pumping although it should drop to zero. The use of non-Markovian master 
equations [2] [17] might overcome these problems, which should also be relevant if we 
ask for the coherence properties of the emitted light [Hill]. Another possibility is the 
combination of CFETs with polynomial techniques for the numerical representation of 
open quantum systems [29H3T] . 

We hope to be able to return to some of these issues soon, which we had to leave 
unresolved here. Presently, we can conclude that CFETs are one promising contribution 
to numerical time propagation of complex non-equilibrium quantum systems and 
warrant further exploration. 
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